Geocodificando y capas base
Analisis espacial con R
Bourdin Elian
19/7/2021
library (tmaptools)
library (RUMBA)
library (tidyverse)
library (sf)
library (ggmap)
library (osmdata)1 Geocodificando
Pero… qué es geocodificar?
Es el proceso de asignar coordenadas (ergo, localización espacial) a uno o varios atributos alfanuméricos que constituyen una dirección. Se requiere de 4 argumentos fundamentales para ello: Calle, Altura o número, Localidad y Estado o Provincia. Otros atributos mejoran la calidad y la precisión: Barrio, País, Código Postal, etc.
Para geocodificar hacemos uso de APIs y servicios geográficos, en este caso, veremos cómo usar el de OpenStreetMap o también conocido como OSM (Nominatim)
1.1 OpenStreetMap
OpenStreetMap es un proyecto colaborativo para crear mapas editables y libres. Los mapas se crean utilizando información geográfica capturada con dispositivos GPS móviles, ortofotografías y otras fuentes libres. Esta cartografía, tanto las imágenes creadas como los datos vectoriales almacenados en su base de datos, se distribuye bajo licencia abierta Licencia Abierta de Bases de Datos.
Los contribuidores más entusiastas mapean barrios completos utilizando herramientas GPS para enviar información local completa, actualizada y precisa a OpenStreetMap. Varias empresas y entidades públicas que producen información geográfica también contribuyen al permitir que sus datos sean incluidos. Existen equipos profesionales de contribuidores que se coordinan para agregar y mantener actualizada información georeferenciada de límites políticos, calles, edificios, negocios y otros puntos de interés.
Toda la información disponible en OpenStreetMap puede ser descargada y reutilizada por cualquier persona, ya sea accediendo al mapa online, obteniendo una copia completa de la base de datos, o accediendo a los datos vía API.
Empecemos!
Vamos a trabajar con datos de una base de datos de entros de salud municipales de Cordoba capital.
setwd("C:/Users/Elian/OneDrive/Escritorio/Ciencia de datos/Geo en Salud/Clase 4") #definiendo el directorio de trabajoCentros_Salud <- readxl::read_xlsx("Listado_Centros_de_Salud_Municipales.xlsx")#deben colocar el directorio en el cual descomprimieron el archivo
#vamos a explorar el objeto
head (Centros_Salud)## # A tibble: 6 x 8
## Z Nº `Centro de Salud` Calle Altura Localidad Horarios Telefono
## <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1 1 GENERAL MOSCONI Pedro Naon 1330 Cordoba 7-21 hs 4335102 /~
## 2 3 2 LOS SAUCES Algarrobo 351 Cordoba 7-19 hs 4338516
## 3 2 3 LOS PINOS Pje. Gabl~ 3786 Cordoba 7-14hs 4339101
## 4 1 4 NUEVA ITALIA Jose de Q~ 2430 Cordoba 7-17 hs 4335123
## 5 2 5 YOFRE NORTE Rene de C~ 1875 Cordoba 7-14 hs 4339103
## 6 6 6 VILLA RIVERA INDA~ Chamico e~ <NA> Cordoba 7-16.12~ 03543-449~summary(Centros_Salud)## Z Nº Centro de Salud Calle
## Min. :1.00 Min. : 1.00 Length:100 Length:100
## 1st Qu.:2.00 1st Qu.: 25.75 Class :character Class :character
## Median :4.00 Median : 51.50 Mode :character Mode :character
## Mean :3.64 Mean : 51.15
## 3rd Qu.:5.00 3rd Qu.: 76.25
## Max. :6.00 Max. :101.00
## Altura Localidad Horarios Telefono
## Length:100 Length:100 Length:100 Length:100
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
## #cuántos valores nulos tenemos en los datos de dirección?
sum(is.na(Centros_Salud$Calle))## [1] 0sum(is.na(Centros_Salud$Altura)) #Tenemos varios valores nulos## [1] 26Para acceder a las funcionalidades de la API usaremos la libería tmaptools, que tiene una amplia funcionalidad, entre ellas, conectarse a los servicios de OSM.
# primero tenemos que construir el campo de dirección tal como lo requerimos
Centros_Salud_Direccion <- Centros_Salud %>%
mutate(DirCompleta = paste(Calle, Altura, ",", Localidad, ", Argentina"))
#vamos a explorar el objeto
head(Centros_Salud_Direccion)## # A tibble: 6 x 9
## Z Nº `Centro de Salu~ Calle Altura Localidad Horarios Telefono
## <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1 1 GENERAL MOSCONI Pedr~ 1330 Cordoba 7-21 hs 4335102~
## 2 3 2 LOS SAUCES Alga~ 351 Cordoba 7-19 hs 4338516
## 3 2 3 LOS PINOS Pje.~ 3786 Cordoba 7-14hs 4339101
## 4 1 4 NUEVA ITALIA Jose~ 2430 Cordoba 7-17 hs 4335123
## 5 2 5 YOFRE NORTE Rene~ 1875 Cordoba 7-14 hs 4339103
## 6 6 6 VILLA RIVERA IN~ Cham~ <NA> Cordoba 7-16.12~ 03543-4~
## # ... with 1 more variable: DirCompleta <chr>summary(Centros_Salud_Direccion)## Z Nº Centro de Salud Calle
## Min. :1.00 Min. : 1.00 Length:100 Length:100
## 1st Qu.:2.00 1st Qu.: 25.75 Class :character Class :character
## Median :4.00 Median : 51.50 Mode :character Mode :character
## Mean :3.64 Mean : 51.15
## 3rd Qu.:5.00 3rd Qu.: 76.25
## Max. :6.00 Max. :101.00
## Altura Localidad Horarios Telefono
## Length:100 Length:100 Length:100 Length:100
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## DirCompleta
## Length:100
## Class :character
## Mode :character
##
##
## #cuántas direcciones tenemos?
nrow (Centros_Salud_Direccion)## [1] 100Vamos a geocodificar nuestras direcciones
Centros_Salud_Direccion <- Centros_Salud_Direccion %>%
mutate(geo=geocode_OSM( #me devuelve las coordenadas (geocode_OSM) en un campo adicional a los originales
Centros_Salud_Direccion$DirCompleta,
return.first.only = TRUE, # me devuelve solo el primer resultado encontrado
keep.unfound = TRUE, # que mantenga aquellas direcciones que No pudo geolocalizar, con un null en el campo "geo"
details = FALSE, # que no devuelva detalles, id de OSM u objetos cercanos
#as.data.frame = NA,
as.sf = FALSE, # que lo deje como una tabla y no lo convierta en un spatial frame
geometry = "point", # que la geometría sea de tipo PUNTO (es decir, un par de coordenadas)
server = "https://nominatim.openstreetmap.org" # el servicio que deseo usar: Nominatim
))
#qué obtuvimos?
head(Centros_Salud_Direccion)## # A tibble: 6 x 10
## Z Nº `Centro de Salu~ Calle Altura Localidad Horarios Telefono
## <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 1 1 GENERAL MOSCONI Pedr~ 1330 Cordoba 7-21 hs 4335102~
## 2 3 2 LOS SAUCES Alga~ 351 Cordoba 7-19 hs 4338516
## 3 2 3 LOS PINOS Pje.~ 3786 Cordoba 7-14hs 4339101
## 4 1 4 NUEVA ITALIA Jose~ 2430 Cordoba 7-17 hs 4335123
## 5 2 5 YOFRE NORTE Rene~ 1875 Cordoba 7-14 hs 4339103
## 6 6 6 VILLA RIVERA IN~ Cham~ <NA> Cordoba 7-16.12~ 03543-4~
## # ... with 8 more variables: DirCompleta <chr>, geo$query <chr>, $lat <dbl>,
## # $lon <dbl>, $lat_min <dbl>, $lat_max <dbl>, $lon_min <dbl>, $lon_max <dbl>summary (Centros_Salud_Direccion$geo) #vamos a explorar el componente geo del objeto que nos devuelve OSM## query lat lon lat_min
## Length:100 Min. :-33.27 Min. :-64.47 Min. :-33.27
## Class :character 1st Qu.:-31.45 1st Qu.:-64.24 1st Qu.:-31.45
## Mode :character Median :-31.42 Median :-64.21 Median :-31.42
## Mean :-31.52 Mean :-64.13 Mean :-31.52
## 3rd Qu.:-31.38 3rd Qu.:-64.14 3rd Qu.:-31.38
## Max. :-31.32 Max. :-62.19 Max. :-31.32
## NA's :55 NA's :55 NA's :55
## lat_max lon_min lon_max
## Min. :-33.27 Min. :-64.47 Min. :-64.47
## 1st Qu.:-31.45 1st Qu.:-64.24 1st Qu.:-64.24
## Median :-31.42 Median :-64.21 Median :-64.20
## Mean :-31.52 Mean :-64.13 Mean :-64.13
## 3rd Qu.:-31.38 3rd Qu.:-64.15 3rd Qu.:-64.14
## Max. :-31.32 Max. :-62.19 Max. :-62.19
## NA's :55 NA's :55 NA's :55# Base final con coordenadas
Centros_Salud_Direccion <- Centros_Salud_Direccion %>%
mutate(lon=geo$lon, lat = geo$lat) %>% # el campo geo se arma como lista, tenemos que seleccionar una lon y una lat
select(`Centro de Salud`, DirCompleta, lon, lat)
#exploremos el objeto final...
summary (Centros_Salud_Direccion)## Centro de Salud DirCompleta lon lat
## Length:100 Length:100 Min. :-64.47 Min. :-33.27
## Class :character Class :character 1st Qu.:-64.24 1st Qu.:-31.45
## Mode :character Mode :character Median :-64.21 Median :-31.42
## Mean :-64.13 Mean :-31.52
## 3rd Qu.:-64.14 3rd Qu.:-31.38
## Max. :-62.19 Max. :-31.32
## NA's :55 NA's :55#qué nos falta para que sea un objeto geográfico? Agregar el sistema de referencia!
Centros_Salud_Direccion <- Centros_Salud_Direccion %>%
filter(!is.na(lon), !is.na(lat)) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326)
summary (Centros_Salud_Direccion)## Centro de Salud DirCompleta geometry
## Length:45 Length:45 POINT :45
## Class :character Class :character epsg:4326 : 0
## Mode :character Mode :character +proj=long...: 0class(Centros_Salud_Direccion)## [1] "sf" "tbl_df" "tbl" "data.frame"1.2 USIG
Hay otros servicios para trabajar datos!
USIG es la Unidad de Sistemas de Información Geográfica del Gobierno de la Ciudad de Buenos Aires y brinda un servicio de geocodificación (entre otras cosas) de datos de la Ciudad de Buenos Aires de manera gratuita!!!
La librería que nos facilita el uso de este servicio se llama RUMBA, la cual la instalamos con el clásico install.packages(“RUMBA”).
RUMBA es un conjunto de herramientas para el análisis de la Región Urbana Metropolitana de Buenos Aires usando R. Incluye funciones que permiten obtener coordenadas precisas (longitud y latitud) que corresponden a direcciones dentro de los límites de la Región Urbana Metropolitana de Buenos Aires o bien de barrios vulnerables de la Ciudad de Buenos Aires.
Las funciones consultan la API del Normalizador de direcciones y de Déficit Habitacional de la USIG. Ademas de las coordenadas, se obtiene la dirección normalizada (escrita de forma inequívoca).
Les proponemos hacer un ejercicio con las postas de vacunación de CABA
Postas_Vacunacion <- readxl::read_xlsx("postas_vacunacion_covid.xlsx")#deben colocar el directorio en el cual descomprimieron el archivo
#vamos a explorar el objeto
head(Postas_Vacunacion)## # A tibble: 6 x 7
## categoria clasif_int efector tipo_de_ef direccion barrio comuna
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 Público Atención a ~ Oficina Na~ Sistema Público~ Libertad 1~ "RETIR~ Comun~
## 2 Público Atención a ~ Sede Oscoe~ Sistema Público~ Bernardo d~ "BERNA~ MONSE~
## 3 Público Atención a ~ Club San L~ Sistema Público~ Av. La Pla~ "BOEDO" Comun~
## 4 Público Atención a ~ Club Itali~ Sistema Público~ Yerbal 150 "CABAL~ Comun~
## 5 Público Atención a ~ Parque Roca Sistema Público~ Av. Cnel. ~ "CNEL.~ VILLA~
## 6 Público Atención a ~ Club Glori~ Sistema Público~ Bragado 68~ "CABA\~ MATAD~summary(Postas_Vacunacion)## categoria clasif_int efector tipo_de_ef
## Length:85 Length:85 Length:85 Length:85
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## direccion barrio comuna
## Length:85 Length:85 Length:85
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :characternrow(Postas_Vacunacion)## [1] 85El dataframe tiene una columna que hace referencia a la dirección de los centros de vacunación. Pero no tenemos ni latitud ni longitud. Ya sabemos que podemos geocodificar los datos! Vamos a probar el servicio de la USIG.
Con la función mutate_USIG_geocode podemos agregar las columnas de lon y lat y así tener las coordenadas. La dirección debe estar expresada como “calle altura, partido”, “calle altura, municipio”, “calle y calle, partido”, o “calle altura, municipio”. El partido o municipio son opcionales. De no ser aclarados, y encontrarse múltiples direcciones que coincidan con la búsqueda, se entregaran las coordenadas dentro de la Ciudad Autónoma de Buenos Aires (si existieran), o en su defecto las del primer partido -por orden alfabético- donde se haya encontrado la dirección.
En resumen: es mejor incluir partido o municipio en las direcciones a georeferenciar.
Postas_Vacunacion_Direccion <- mutate_USIG_geocode(Postas_Vacunacion, "direccion") #paciencia. Puede tardar un poco! Utilizar la función es tan simple como indicar el DF y el nombre del campo que tiene la dirección que vamos a trabajar.
#vamos a explorar el objeto
head(Postas_Vacunacion_Direccion)## categoria clasif_int efector
## 1 Público Atención a Adultos Mayores Oficina Nacional Scout
## 2 Público Atención a Adultos Mayores Sede Oscoema
## 3 Público Atención a Adultos Mayores Club San Lorenzo (Sede Av. La Plata)
## 4 Público Atención a Adultos Mayores Club Italiano
## 5 Público Atención a Adultos Mayores Parque Roca
## 6 Público Atención a Adultos Mayores Club Glorias Argentinas
## tipo_de_ef
## 1 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 2 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 3 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 4 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 5 Sistema Público de la Ciudad (Posta extrahospitalaria)
## 6 Sistema Público de la Ciudad (Posta extrahospitalaria)
## direccion barrio comuna
## 1 Libertad 1282 RETIRO Comuna 1
## 2 Bernardo de Irigoyen 330 BERNARDO DE 330" MONSERRAT
## 3 Av. La Plata 1770 BOEDO Comuna 5
## 4 Yerbal 150 CABALLITO Comuna 6
## 5 Av. Cnel. Roca 4200 CNEL. AV. 4200" VILLA SOLDATI
## 6 Bragado 6875 CABA" MATADEROS
## address_normalised lon lat
## 1 LIBERTAD 1282, CABA -58.38452 -34.59320
## 2 IRIGOYEN, BERNARDO DE 330, CABA -58.38049 -34.61215
## 3 LA PLATA AV. 1770, CABA -58.42459 -34.63538
## 4 YERBAL 150, CABA -58.43292 -34.61538
## 5 ROCA, CNEL. AV. 4200, CABA -58.44747 -34.67549
## 6 BRAGADO 6875, CABA -58.50989 -34.65928summary(Postas_Vacunacion_Direccion)## categoria clasif_int efector tipo_de_ef
## Length:85 Length:85 Length:85 Length:85
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## direccion barrio comuna address_normalised
## Length:85 Length:85 Length:85 Length:85
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## lon lat
## Min. :-58.52 Min. :-34.68
## 1st Qu.:-58.47 1st Qu.:-34.63
## Median :-58.43 Median :-34.62
## Mean :-58.43 Mean :-34.62
## 3rd Qu.:-58.40 3rd Qu.:-34.59
## Max. :-58.36 Max. :-34.55
## NA's :33 NA's :33#Solo nos queda transformarlo a objeto geográfico
Postas_Vacunacion_Direccion <- Postas_Vacunacion_Direccion %>%
filter(!is.na(lon), !is.na(lat)) %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326)
summary (Postas_Vacunacion_Direccion)## categoria clasif_int efector tipo_de_ef
## Length:52 Length:52 Length:52 Length:52
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## direccion barrio comuna address_normalised
## Length:52 Length:52 Length:52 Length:52
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## geometry
## POINT :52
## epsg:4326 : 0
## +proj=long...: 0class(Postas_Vacunacion_Direccion)## [1] "sf" "data.frame"2 Agregando contexto a la información
Pudimos geolocalizar gran parte de nuestras direcciones de interés, aunque ya vimos que no es un proceso infalible… pero cuando trabajamos con datos gegráficos es importante poder tener una referencia, ubicar los puntos en contexto, pero ¿Qué pasa cuando no tenemos un mapa?
Tenemos una opción! Podemos agregar capas base, que sería como agregar “una foto” sobre la cual posicionaremos los datos de interés.
El primer paso para obtener las capas bases que necesitamos es obtener el bounding box del mapa que queremos descargar. ¿Qué utilizamos?… así es! getbb. En este caso obtendremos un mapa de la Ciudad Autónoma de Buenos Aires, ya que vamos a localizar las Postas de Vacunación que geocodificamos en el paso anterior.
bbox <- getbb("Ciudad Autónoma de Buenos Aires,Argentina")
head(bbox)## min max
## x -58.53145 -58.33514
## y -34.70564 -34.52655Ya tenemos los límites de la ciudad en la variable bbox. Ahora vamos a descargar la capa base utilizando ggmap. Lo que permite hacer ggmap es, en esencia, es añadir a los gráficos ya conocidos una capa cartográfica adicional. Para eso usa recursos disponibles en la web a través de APIs; ggmap tiene funciones para obtener mapas (de diversos tipos y de distintos orígenes: Google, Stamen, OpenStreetMap). Nosotros vamos a obtener un mapa de Stamen que es un estudio de visualización de datos que abre sus mapas para el uso de la comunidad. Pueden ver sus diferentes diseños en Stamen maps
#Vamos a descargar el mapa
CABATerra <- get_stamenmap(bbox = bbox, #En bboxle indicamos los límites específicos, ya lo sabemos!
maptype = "terrain", #En maptype le indicamos el tipo de mapa que queremos descargarnos.
zoom=13)#En zoom le indicamos el nivel de detalle que queremos. A más zoom, más detalle, pero también más pesado el archivo que nos descarguemos…
#que tenemos en la variable CABA?
head(CABATerra) #raro?## [1] "#EAEAEA" "#FAFAFA" "#D1D1D1" "#999999" "#999999" "#999999"CABATerra #R ya nos dice en la descripción que estamos trabajando con una imagen de mapa de Stamen Maps... e incluso nos sugiere como graficarlo!!!## 1268x1143 terrain map image from Stamen Maps.
## See ?ggmap to plot it.Vamos a hacerle caso a los consejos de R y vamos a utilizar ggmap para dibujar el mapa. Por suerte para nosotros, ggmap tiene la misma sintaxis que ggplot(), por lo tanto no tendremos demasiados inconvenientes.
ggmap(CABATerra)
Y si probamos con otro tipo de mapa?
CABATonner <- get_stamenmap(bbox = bbox,
maptype = "toner-background",
zoom=13)
ggmap(CABATonner)
Excelente, tenemos nuestra capa de base, pero ahora cómo ubicamos en contexto las Postas de Vacunación de la Ciudad? Fácil… agregando una nueva capa a nuestro mapa!
ggmap(CABATonner) +
geom_sf(data=Postas_Vacunacion_Direccion, inherit.aes=FALSE, color = "red", size = 2)## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
¿Qué les parece el mapa que obtuvimos?…